#include "mesh.hpp"
#include "util.hpp"
#include "meshFunctions.hpp"
#include "mesh_testFunctions.hpp"
#include "testFunctions.hpp"

#include <catch2/catch_test_macros.hpp>

using namespace hms;

TEST_CASE( "Mesh::findDomainCellEdges" )
{
	SECTION( "rectangular mesh" )
	{
		auto mesh { testMeshRectangles() };
		auto domainCellEdges { findDomainCellEdges(mesh) };
		
		REQUIRE( sameEntries(
			domainCellEdges.second.col(0),
			IdxVector { 0, 2, 6, 8 }
		));
		REQUIRE( sameEntries(
			domainCellEdges.second.col(1),
			IdxVector { 1, 3, 8,10 }
		));
		REQUIRE( sameEntries(
			domainCellEdges.second.col(2),
			IdxVector { 2, 4, 7, 9 }
		));
		REQUIRE( sameEntries(
			domainCellEdges.second.col(3),
			IdxVector { 3, 5, 9,11 }
		));
	}
	SECTION( "triangular mesh" )
	{
		auto mesh { testMeshTriangles() };
		auto domainCellEdges { findDomainCellEdges(mesh) };
		
		REQUIRE( sameEntries(
			domainCellEdges.second.col(0).head(3),
			IdxVector { 0, 1, 8 }
		));
		REQUIRE( sameEntries(
			domainCellEdges.second.col(1).head(3),
			IdxVector { 1, 2, 7 }
		));
		REQUIRE( sameEntries(
			domainCellEdges.second.col(2).head(3),
			IdxVector { 2, 3, 6 }
		));
		REQUIRE( sameEntries(
			domainCellEdges.second.col(3).head(3),
			IdxVector { 3, 4, 5 }
		));
	}
	SECTION( "mixed mesh" )
	{
		auto mesh { testMeshMixed() };
		auto domainCellEdges { findDomainCellEdges(mesh) };
		
		REQUIRE( sameEntries(
			domainCellEdges.second.col(0).head(3),
			IdxVector { 0, 1, 2 }
		));
		REQUIRE( sameEntries(
			domainCellEdges.second.col(1),
			IdxVector { 2, 3, 4, 5 }
		));
	}
}

 
TEST_CASE( "Swapping edges" ){
	auto mesh { testMeshRectangles(true) };
	Array2Xi
		edgeVertexIndices { mesh.edges().vertices() }, vertsCopy,
		edgeAdjacentCells { mesh.edges().cells   () }, cellsCopy;
	Array4Xi domainCellEdges { mesh.cells().domainCellEdges() }, edgesCopy;

	vertsCopy = edgeVertexIndices;
	cellsCopy = edgeAdjacentCells;
	edgesCopy = domainCellEdges;

	SECTION( "Swapping only vertex and parent cell indices" ){
		swapEdges(0, 8, edgeVertexIndices, edgeAdjacentCells);

		REQUIRE( (edgeVertexIndices.col(0) == vertsCopy.col(8)).all() );
		REQUIRE( (edgeVertexIndices.col(8) == vertsCopy.col(0)).all() );
		REQUIRE( (edgeAdjacentCells.col(0) == cellsCopy.col(8)).all() );
		REQUIRE( (edgeAdjacentCells.col(8) == cellsCopy.col(0)).all() );
	}
	SECTION( "Swapping vertex and parent cell indices "
		"as well as domain cell edge indices"
	){
		swapEdges(0, 8, edgeVertexIndices, edgeAdjacentCells, domainCellEdges);

		REQUIRE( (edgeVertexIndices.col(0) == vertsCopy.col(8)).all() );
		REQUIRE( (edgeVertexIndices.col(8) == vertsCopy.col(0)).all() );
		REQUIRE( (edgeAdjacentCells.col(0) == cellsCopy.col(8)).all() );
		REQUIRE( (edgeAdjacentCells.col(8) == cellsCopy.col(0)).all() );
		REQUIRE( (domainCellEdges.col(0) == Array4i{0, 8, 2, 6}).all() );
		REQUIRE( (domainCellEdges.col(1) == Array4i{9,10, 3, 8}).all() );
		REQUIRE( (domainCellEdges.rightCols(2) == edgesCopy.rightCols(2)).all() );
		
		swapEdges(1, 9, edgeVertexIndices, edgeAdjacentCells, domainCellEdges);
		REQUIRE( (domainCellEdges.col(0) == Array4i{0, 8, 2, 6}).all() );
		REQUIRE( (domainCellEdges.col(1) == Array4i{1,10, 3, 8}).all() );
		REQUIRE( (domainCellEdges.col(2) == Array4i{2, 9, 4, 7}).all() );
		REQUIRE( (domainCellEdges.col(3) == Array4i{3,11, 5, 9}).all() );
		
		swapEdges(2, 3, edgeVertexIndices, edgeAdjacentCells, domainCellEdges);
		REQUIRE( (domainCellEdges.col(0) == Array4i{0, 8, 3, 6}).all() );
		REQUIRE( (domainCellEdges.col(1) == Array4i{1,10, 2, 8}).all() );
		REQUIRE( (domainCellEdges.col(2) == Array4i{3, 9, 4, 7}).all() );
		REQUIRE( (domainCellEdges.col(3) == Array4i{2,11, 5, 9}).all() );
	}

	// std::cout
	// 	<<   "edgeVertexIndices (before):\n" << vertsCopy
	// 	<< "\nedgeVertexIndices (after): \n" << edgeVertexIndices
	// 	<< "\nedgeAdjacentCells (before):\n" << cellsCopy
	// 	<< "\nedgeAdjacentCells (after): \n" << edgeAdjacentCells
	// 	<< "\ndomainCellEdges   (before):\n" << edgesCopy
	// 	<< "\ndomainCellEdges   (after): \n" << domainCellEdges
	// 	<< "\n";
}

/* making finding boundary edges a bit more convenient */
IdxVector getBoundaryEdges( const hms::Mesh& mesh ){
	std::pair domainCells { hms::findDomainCellEdges(mesh) };
	
	Array2Xi adjacentCells { hms::findAdjacentCells(
		domainCells.first,
		domainCells.second,
		mesh.nEdges()
	)};
	
	return hms::findEdgesWEmptyParent(adjacentCells);
}
	

TEST_CASE( "Finding boundary edges" )
{
	SECTION( "rectangular mesh" ){
		hms::Mesh mesh { testMeshRectangles() };
		
		IdxVector boundaryEdges { getBoundaryEdges(mesh) };
		
		REQUIRE( sameEntries( boundaryEdges, { 0, 1, 4, 5, 6, 7, 10, 11 } ) );
	}
	SECTION( "triangular mesh" ){
		hms::Mesh mesh { testMeshTriangles() };
		
		IdxVector boundaryEdges { getBoundaryEdges(mesh) };
		REQUIRE( sameEntries( boundaryEdges, { 0, 4, 5, 6, 7, 8 } ) );
	}
	SECTION( "mixed mesh" ){
		hms::Mesh mesh { testMeshMixed() };
		
		IdxVector boundaryEdges { getBoundaryEdges(mesh) };
		REQUIRE( sameEntries( boundaryEdges, { 0, 1, 3, 4, 5 } ) );
	}
}

TEST_CASE( "Creating cells from edges" )
{
	SECTION( "rectangular mesh" ) {
		hms::Mesh mesh { testMeshRectangles(true) };

		const auto& cellEdges  = mesh.cells().domainCellEdges();
		const auto& nEdges = mesh.cells().nDomainCellEdges();

		REQUIRE( sameEntries(cellEdges.col(0), Array4i{0, 2, 6, 8}) );
		REQUIRE( sameEntries(cellEdges.col(1), Array4i{0, 3, 9,10}) );
		REQUIRE( sameEntries(cellEdges.col(2), Array4i{1, 2, 4, 7}) );
		REQUIRE( sameEntries(cellEdges.col(3), Array4i{1, 3, 5,11}) );
		
		REQUIRE( (nEdges == Array4<char>{4,4,4,4}).all() );
	}
	SECTION( "triangular mesh" ) {
		hms::Mesh mesh { testMeshTriangles(true) };

		const auto& cellEdges  = mesh.cells().domainCellEdges();
		const auto& nEdges = mesh.cells().nDomainCellEdges();

		REQUIRE( sameEntries(
			cellEdges.col(0).head(3), Vector3i{1, 3, 8} ) );
		REQUIRE( sameEntries(
			cellEdges.col(1).head(3), Vector3i{1, 2, 7} ) );
		REQUIRE( sameEntries(
			cellEdges.col(2).head(3), Vector3i{0, 2, 6} ) );
		REQUIRE( sameEntries(
			cellEdges.col(3).head(3), Vector3i{0, 4, 5} ) );
		
		REQUIRE( (nEdges == Array4<char>{3,3,3,3}).all() );
	}
	SECTION( "mixed mesh" ) {
		hms::Mesh mesh { testMeshMixed(true) };

		const auto& cellEdges  = mesh.cells().domainCellEdges();
		const auto& nEdges = mesh.cells().nDomainCellEdges();

		REQUIRE( sameEntries(
			cellEdges.col(0).head(3), Vector3i{0, 1, 2} ) );
		REQUIRE( sameEntries(cellEdges.col(1), Array4i{0, 3, 4, 5}) );
		
		REQUIRE( (nEdges == Array2<char>{3,4}).all() );
	}
}

#include "topology.hpp"
TEST_CASE( "Creating rectangular mesh" ){
	
	auto mesh { hms::createRectangularMesh( 2,2,1,1,-1,-1 ) };
	
	REQUIRE( mesh.vertex(0) == Vector2s{-1,-1} );
	REQUIRE( mesh.vertex(1) == Vector2s{ 0,-1} );
	REQUIRE( mesh.vertex(2) == Vector2s{ 1,-1} );
	REQUIRE( mesh.vertex(3) == Vector2s{-1, 0} );
	REQUIRE( mesh.vertex(4) == Vector2s{ 0, 0} );
	REQUIRE( mesh.vertex(5) == Vector2s{ 1, 0} );
	REQUIRE( mesh.vertex(6) == Vector2s{-1, 1} );
	REQUIRE( mesh.vertex(7) == Vector2s{ 0, 1} );
	REQUIRE( mesh.vertex(8) == Vector2s{ 1, 1} );
	
	/* vertices must be ordered counter-clockwise */
	REQUIRE( (mesh.cell(0).vertices(true) == Array4i {0,1,4,3}).all() );
	REQUIRE( (mesh.cell(1).vertices(true) == Array4i {1,2,5,4}).all() );
	REQUIRE( (mesh.cell(2).vertices(true) == Array4i {3,4,7,6}).all() );
	REQUIRE( (mesh.cell(3).vertices(true) == Array4i {4,5,8,7}).all() );
	
	/* the edge indexes aren't actually a requirement */
	// REQUIRE( sameEntries(mesh.cells().edges(0), IdxVector { 1, 2, 6, 7}) );
	// REQUIRE( sameEntries(mesh.cells().edges(1), IdxVector { 0, 1, 4, 5}) );
	// REQUIRE( sameEntries(mesh.cells().edges(2), IdxVector { 2, 3, 9, 8}) );
	// REQUIRE( sameEntries(mesh.cells().edges(3), IdxVector { 0, 3,10,11}) );
}


TEST_CASE( "Mesh from vertices and cells", "[only]" ){
	Index nDomainCells {6};
	Array2Xs vertices {
		ArrayNs<3*12>{
			0, 0, 0, 1.1, 0, 0, 2.2, 0, 0, 0, 1.4, 0, 1.1, 1.4, 0, 2.2, 1.4, 0,
			0, 2.8, 0, 1.1, 2.8, 0, 2.2, 2.8, 0, 0, 4.2, 0, 1.1, 4.2, 0, 2.2, 4.2, 0
		}.reshaped(3,12).topRows<2>()
	};
	Array4Xi domainCellVertices {
		ArrayNi<24>{
			0, 1, 4, 3, 1, 2, 5, 4, 3, 4, 7, 6,
			4, 5, 8, 7, 6, 7, 10, 9, 7, 8, 11, 10
		}.reshaped(4, nDomainCells)
	};
	ArrayX<char> nDomainCellEdges (nDomainCells);
	nDomainCellEdges = 4;

	auto mesh { cellVerticesToMesh(
		std::move(vertices), std::move(nDomainCellEdges), domainCellVertices
	) };

	// Array4Xi domainCellVerticesStored ( 4, mesh.nDomainCells() );
	// ArrayNXs<8> domainCellVertexCoords ( 8, mesh.nDomainCells() );
	// for ( int i=0; i<mesh.nDomainCells(); ++i ){
	// 	domainCellVerticesStored.col(i) = mesh.cell(i).vertices(true);
	// 	domainCellVertexCoords.col(i) = mesh.vertices()(
	// 		hms::all, domainCellVerticesStored.col(i)
	// 	).reshaped();
	// }

	// std::cout << "Vertices:\n" << mesh.vertices()
	// 	<< "\nEdge vertex indices:\n" << mesh.edges().vertices()
	// 	<< "\nEdge lengths:\n" << mesh.edges().lengths().transpose()
	// 	<< "\nEdge adjacent cells:\n" << mesh.edges().cells()
	// 	<< "\nCell vertex indices:\n" << domainCellVerticesStored
	// 	<< "\nCell vertex coords:\n" << domainCellVertexCoords
	// 	<< "\nCell edge indices:\n" << mesh.cells().domainCellEdges()
	// 	<< "\nCell areas:\n" << mesh.cells().areas().transpose()
	// 	<< "\nCell centroids:\n" << mesh.cells().centroids()
	// 	<< "\nCell circumfs:\n" << mesh.cells().circumfs().transpose()
	// 	<< "\n";
}
